#### REPLICATION MATERIALS FOR ''BENCHMARKING PANDEMIC RESPONSE: HOW THE UK'S COVID-19 VACCINE ROLL-OUT AFFECTED POPULAR SUPPORT FOR THE EU''
###### Authors: Irene Rodríguez, Toni Rodon, Asli Unan, Lisa Herbig, Heike Klüver & Theresa Kuhn
###### Journal: British Journal of Political Science

###### FIGURES AND TABLES IN APPENDICES

# Loading libraries
if (!require("pacman")) install.packages("pacman")
pacman::p_load(kableExtra, marginaleffects, modelsummary, knitr, estimatr, tidyverse, cobalt, knitr, 
               RColorBrewer, janitor, lubridate, readstata13, countrycode, rdrobust, magick, jtools, dplyr,
               ggiraphExtra, sjPlot, gridExtra, data.table, pastecs, MatchIt, psych, gridExtra, webshot2)

# Set working directory

setwd("/Replication materials")

load("Data/data.RData")


#### Table 1: Dependent variables descriptives ####

# Selecting dependent variables
data_dep_1 <- data |> 
  select(eu_support, eu_positive, eu_benefits, eu_right_direction, eu_democracy, eu_help, eu_coord, eu_health_issues,
         eu_ss_issues, eu_education_issues, eu_equality_issues, eu_climate_issues, eu_jobs_issues, eu_digitalisation_issues, eu_workingconditions_issues) 

# Extracting descriptive measures
descriptives_dependent <- t(data.frame(stat.desc(data_dep_1))) |> 
  as.data.frame() |> 
  rownames_to_column() |> 
  select(!"nbr.null")

# Table containing descriptive measures of dependent variables
t_descriptives_dependent <- kbl(descriptives_dependent, 
                                row.names = FALSE,
                                escape = FALSE,
                                digits = 2)

# Saving Table 2 as app_table_1.png
save_kable(
  t_descriptives_dependent,
  "Tables/app_table_1.png",
  bs_theme = "journal",
  density = 600
)

#### Table 2: Numeric independent variables descriptives ####

# Selecting numerical independent variables
dat_ind_num <- data |> 
  select(age, eduy, lr_self, local_interest, national_interest, european_interest, hospital_occupancy, mass_ideo, elite_ideo)

# Extracting descriptive measures
descriptives_independent <- t(data.frame(stat.desc(dat_ind_num))) |> 
  as.data.frame() |> 
  rownames_to_column() |> 
  select(!"nbr.null")

# Table containing descriptive masures of dependent variables

t_descriptives_independent <- kbl(descriptives_independent, 
                                  row.names = FALSE,
                                  escape = FALSE,
                                  digits = 2)

# Saving Table 2 as app_table_2.png
save_kable(
  t_descriptives_independent,
  "~/Replication materials/Tables/app_table_2.png",
  bs_theme = "journal",
  density = 600)

#### Figure 1: Categorical independent variables descriptives  ####

# Selecting independent categorical variables from the dataframe
dat_ind_cat <- data |> 
  select(employment, living_area, gender, ideology_r, treat2, europe)

### plot for employment
plot_employ <- ggplot(data = subset(dat_ind_cat, 
                                    !is.na(employment)), 
                      mapping = aes(x = employment)) +
  geom_bar() +
  theme_classic() +
  labs(y = "Frequency",
       x = "Employment") 

### plot for living area
plot_rural <- ggplot(data = subset(dat_ind_cat, !is.na(living_area)), mapping = aes(x = living_area)) +
  geom_bar() +
  scale_x_discrete(limits = c(1, 2, 3),
                   breaks = c(1, 2, 3),
                   labels = c("Rural area", "Small/middle town", "Large town")) +
  theme_classic() +
  
  labs(y = "Frequency",
       x = "Rural / urban") 

### plot for gender
plot_gender <- ggplot(data = subset(dat_ind_cat, !is.na(gender)), mapping = aes(x = gender)) +
  geom_bar() +
  scale_x_discrete(limits = c(0, 1),
                   breaks = c(0, 1),
                   labels = c("Female", "Male")) +
  theme_classic() +
  labs(y = "Frequency",
       x = "Gender") 

### plot for left-right self-placement
plot_ideology <- ggplot(data = subset(dat_ind_cat, 
                                      !is.na(ideology_r)),
                        mapping = aes(x = ideology_r)) +
  geom_bar() +
  theme_classic() +
  scale_x_discrete(limits = c("far left", 
                              "center-left", 
                              "center", 
                              "center_right", 
                              "far right"),
                   breaks = c("far left", 
                              "center-left", 
                              "center", 
                              "center_right", 
                              "far right"),
                   labels = c("Far left", 
                              "Center left", 
                              "Center", 
                              "Center right", 
                              "Far right")) +
  labs(y = "Frequency",
       x = "Ideology") 

### plot for treatment
plot_treat <- ggplot(data = subset(dat_ind_cat, 
                                   !is.na(treat2)), 
                     mapping = aes(x = treat2)) +
  geom_bar() +
  scale_x_discrete(limits = c(0, 1),
                   breaks = c(0, 1),
                   labels = c("Before treatment", "After treatment")) +
  theme_classic() +
  labs(y = "Frequency",
       x = "Treatment")

### plot for region
plot_region <- ggplot(data = subset(dat_ind_cat, 
                                    !is.na(europe)), 
                      mapping = aes(x = europe)) +
  geom_bar() +
  theme_classic() +
  scale_x_discrete(labels = c("Central", "Eastern", "Northern", "Southern", "Western")) +
  labs(y = "Frequency",
       x = "Region")

### merging all plots into one figure
plots_descriptives <- grid.arrange(arrangeGrob(plot_employ, plot_gender, plot_rural, ncol = 3),
                                   arrangeGrob(plot_ideology, plot_treat, plot_region, ncol = 3),
                                   nrow = 2)

ggsave("~/Replication materials/Figures/app_figure_1.png", 
     plot = plots_descriptives, 
       dpi = 600, 
       width = 11, 
       height = 5)

#### Table 3: Newspapers by country ####

# Adding data to a data frame
newspapers_by_country <- data.frame(
  Country = c("Germany", "Germany", "Germany", "The Netherlands", "The Netherlands", "The Netherlands","Poland", "Poland", "Poland", "Spain", "Spain", "Spain"),
  Newspaper = c("Bild", "Die Welt", "Süddeutsche Zeitung", "De Telegraaf", "NRC", "Volkskrant", "Fakt", "Gazeta Wyborcza", "Rzeczpospolita", "ABC", "El Mundo", "El País"),
  Count = c(126753, 24490, 364212, 111549, 29457, 37600, 50886, 16380, 11987, 31354, 25453, 38738)
)

# Building the table
t_newspapers_by_country <- kbl(newspapers_by_country, 
                               row.names = FALSE,
                               escape = FALSE)

# Saving Table 3 as app_table_3.png
save_kable(
  t_newspapers_by_country,
  "~/Replication materials/Tables/app_table_3.png",
  bs_theme = "journal",
  density = 600)

#### Figure 2: Attention index to the Brexit issue in press articles from European Union countries (weekly) ####

# Data obtained through python scraping of newspaper titles in Spain, Germany, the Netherlands, and Poland
## Newspaper articles that mention Brexit in the media

## Importing dataset
brexit_att_w <- read_csv("Data/brexit_attention_weekly.csv") 

## Manipulating dataset to fit visualisation needs
brexit_att_w <- brexit_att_w |> 
  pivot_wider(names_from = "topic",
              values_from = "attention")

## Building plot for Figure C.1
p_attention_brexit <- ggplot(data = brexit_att_w,
                             mapping = aes(x = date, y = Brexit)) +
  geom_line(color = "blue") +
  labs(x = 'Date',
       y = 'Issue attention',
       color = '',
       title = '') +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 13),
        legend.position = "bottom")

## Saving Figure 2 as 'app_figure_2.png'
ggsave("~/Replication materials/Figures/app_figure_2.png", 
       plot = p_attention_brexit, 
       dpi = 600, 
       width = 9, 
       height = 5)



#### Figure 3: Worldwide searches of the word "Brexit" in Google ####

# Importing dataset
brexit_google <- read_csv("Data/google_brexit.csv",
                          skip = 1,
                          col_names = c("date", "index"),
                          col_types = c("D", "i"))

# Manipulating dataset to fit visualisation needs
brexit_google <- brexit_google[-c(1), ] |> 
  mutate(index = as.numeric(index))

# Building plot for Figure 3
p_google_brexit <- ggplot(data = brexit_google, 
                          mapping = aes(x = date, y = index)) +
  geom_line(color = "blue") +
  labs(x = 'Date',
       y = 'Google searches index',
       color = '') +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 13),
        legend.position = "bottom")

## Saving Figure 3 as 'app_figure_3.png'
ggsave("~/Replication materials/Figures/app_figure_3.png", 
       plot = p_google_brexit, 
       dpi = 600, 
       width = 9, 
       height = 5)

#### Table 4: Treatment time and fieldwork per country in sample ####

## Building data frame
tt_by_country <- data.frame(
  Country = c("Austria", "Bulgaria", "Cyprus", "France", "Germany", "Ireland", "Italy", "Luxembourg", "The Netherlands", "Poland", "Portugal", "Spain"),
  `Treatment date` = c("7 December 2020", "8 December 2020", "7 December 2020", "7 December 2020", "9 December 2020", "7 December 2020", "8 December 2020", "9 December 2020", "9 December 2020", "8 December 2020", "8 December 2020", "7 December 2020"),
  `Fieldwork start` = c("21 November 2020", "24 November 2020", "21 November 2020", "23 November 2020", "24 November 2020", "03 December 2020", "23 November 2020", "03 December 2020", "23 November 2020", "28 November 2020", "26 November 2020", "24 November 2020"),
  `Fieldwork end` = c("10 December 2020", "20 December 2020", "12 December 2020", "16 December 2020", "21 December 2020", "21 December 2020", "15 December 2020", "21 December 2020", "16 December 2020", "14 December 2020", "21 December 2020", "15 December 2020")
)

## Building table for saving
t_treatment_by_country <- kbl(tt_by_country, 
                              row.names = FALSE,
                              escape = FALSE)

## Saving Table 4 as 'app_table_4.png'
save_kable(
  t_treatment_by_country,
  "~/Replication materials/Tables/app_table_4.png",
  bs_theme = "journal",
  density = 600
)

#### Figure 4: Interviews per day and country before and after the start date of the UK vaccine rollout ####

## Building plot
p_fieldwork_countries <- ggplot(data = filter(data, is.na(data$treat2) == F), aes(x=day)) + 
  geom_bar() +
  scale_x_continuous(limits = c(0, 31),
                     breaks = seq(0,31,5)) +
  geom_vline(xintercept = 17,
             color =  "red",
             size = 0.5,
             linetype = "dashed") +
  facet_wrap(~country) +
  labs(y = "Participants",
       x = "Day") +
  theme_classic() +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 13),
        legend.position = "bottom")

## Saving Figure 4 plot as 'app_figure_4.png'
ggsave("~/Replication materials/Figures/app_figure_4.png", 
       plot = p_fieldwork_countries, 
       dpi = 600, 
       width = 8, 
       height = 7)


#### Figure 5: Balance distributions for covariates depending on treatment ####

## Filter the dataset 'dat_2' to select the relevant variables of interest for analysis and omit any rows with missing values (NA)
data_matching <- data |>  
  select(c("treat2", "country", "age", "gender", "eduy", "lr_self", "subjective_class", "living_area", "local_interest", "national_interest", "european_interest")) %>%
  na.omit()

## Matching observations by specific variables
m.out <- MatchIt::matchit(treat2 ~ age + gender + eduy + lr_self + subjective_class + living_area + local_interest + national_interest + european_interest + factor(country),
                          data = data_matching) # Perform propensity score matching using the 'matchit' function

# Assess balance of covariate 'country' post-matching using the 'bal.tab' function
country <- bal.plot(m.out, var.name = "factor(country)", 
                    which = "unadjusted") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1)) +
  labs(title = "Country of the respondent",
       x = "Country")

# Assess balance of covariate 'eduy' post-matching using the 'bal.tab' function
educ <- bal.plot(m.out, var.name = "eduy",
                 which = "unadjusted") +
  theme(legend.position = "none") +
  labs(title = "Education",
       x = "Years of education")

# Assess balance of covariate 'age' post-matching using the 'bal.tab' function
age <- bal.plot(m.out, var.name = "age",
                mirror = TRUE,
                type = "histogram",
                which = "unadjusted") +
  theme(legend.position = "none") +
  labs(title = "Age",
       x = "Years")

# Assess balance of covariate 'gender' post-matching using the 'bal.tab' function
gender <- bal.plot(m.out, var.name = "gender", 
                   type = "histogram",
                   which = "unadjusted") +
  theme(legend.position = "none") +
  labs(title = "Gender",
       x = "") +
  scale_x_discrete(breaks = c(0,1),
                   labels = c("Female", "Male"))

# Assess balance of covariate 'lr_self' post-matching using the 'bal.tab' function
lr <- bal.plot(m.out, var.name = "lr_self",
               type = "histogram",
               which = "unadjusted",
               mirror = TRUE) +
  theme(legend.position = "none") +
  labs(title = "Left-right self-placement",
       x = "Left-right self-placement")

# Assess balance of covariate 'subjective_class' post-matching using the 'bal.tab' function
class <- bal.plot(m.out, var.name = "subjective_class",
                  which = "unadjusted",
                  type = "histogram",
                  mirror = TRUE) +
  theme(legend.position = "none") +
  labs(title = "Subjective class",
       x = "") +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5),
                     labels = c("Lower", "2", "3", "4", "Higher"))

# Assess balance of covariate 'living_area' post-matching using the 'bal.tab' function
rural <- bal.plot(m.out, var.name = "living_area",
                  which = "unadjusted",
                  type = "histogram",
                  mirror = TRUE) +
  labs(title = "Rural/urban cleavage",
       x = "") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = c(1, 2, 3),
                     labels = c("Rural area", "Small/middle town", "Large town"))

# Assess balance of covariate 'local_interest' post-matching using the 'bal.tab' function
localinterest <- bal.plot(m.out, var.name = "local_interest",
                          which = "unadjusted") +
  labs(title = "Interest in local politics",
       x = "Level of interest") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9),
                     labels = c("Low", 2, 3, 4, 5, 6, 7, 8, "High"))

# Assess balance of covariate 'national_interest' post-matching using the 'bal.tab' function
nationalinterest <- bal.plot(m.out, var.name = "national_interest",
                             which = "unadjusted") +
  labs(title = "Interest in national politics",
       x = "Level of interest") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9),
                     labels = c("Low", 2, 3, 4, 5, 6, 7, 8, "High"))

# Assess balance of covariate 'european_interest' post-matching using the 'bal.tab' function
europeaninterest <- bal.plot(m.out, var.name = "european_interest",
                             which = "unadjusted") +
  labs(title = "Interest in European politics",
       x = "Level of interest") +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7, 8, 9),
                     labels = c("Low", 2, 3, 4, 5, 6, 7, 8, "High"))

## Arrange all plots as a single figure
balance <- grid.arrange(educ, age, 
                        gender, lr, 
                        class, localinterest, 
                        nationalinterest, europeaninterest, 
                        rural, country, 
                        nrow = 5)

## Saving Figure 5 as 'app_figure_5.png'
ggsave("~/Replication materials/Figures/app_figure_5.png", 
       plot = balance, 
       dpi = 600, 
       width = 6, 
       height = 10)

#### Table 5: Effect of the UK vaccine rollout on EU attitudes ####

## Model building for DVs: EU support, EU benefits, EU positive image, EU right direction, EU democracy, EU health issues, EU coordination, and EU future coordination
model_eu_support_fe <- lm_robust(eu_support ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                 data = data,
                                 fixed_effects = country, ## adding country fixed effects
                                 se_type = "stata") ## adding robust standard errors
model_eu_beneficial_fe <- lm_robust(eu_benefits ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                    data = data,
                                    fixed_effects = country,
                                    se_type = "stata")
model_eu_positive_fe <- lm_robust(eu_positive ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                  data = data,
                                  fixed_effects = country,
                                  se_type = "stata")
model_eu_right_direction_fe <- lm_robust(eu_right_direction ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                         data = data,
                                         fixed_effects = country,
                                         se_type = "stata")
model_eu_democracy_fe <- lm_robust(eu_democracy ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                   data = data,
                                   fixed_effects = country,
                                   se_type = "stata")
model_eu_health_issues_fe <- lm_robust(eu_health_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                       data = data,
                                       fixed_effects = country,
                                       se_type = "stata")
model_eu_future_coordination_fe <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                             data = data,
                                             fixed_effects = country,
                                             se_type = "stata")
model_eu_coordination_fe <- lm_robust(eu_help ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                      data = data,
                                      fixed_effects = country,
                                      se_type = "stata")

## Making a list of all main models
models_fe <- list("EU support" = model_eu_support_fe,
                  "EU beneficial" = model_eu_beneficial_fe,
                  "Positive image" = model_eu_positive_fe,
                  "EU right direction" = model_eu_right_direction_fe,
                  "EU democracy" = model_eu_democracy_fe,
                  "Health decisions" = model_eu_health_issues_fe,
                  "EU crisis coordination" = model_eu_coordination_fe,
                  "EU future crisis coordination" = model_eu_future_coordination_fe)

## Regression table of main models and saving Table 6 as 'app_table_5.png'
t_models_main <- modelsummary(models_fe,
                              stars = TRUE, ## adding stars for significance
                              coef_map = c('treat2' = 'Treatment',
                                           'gender' = 'Gender',
                                           'age' = 'Age',
                                           'eduy' = 'Education',
                                           'lr_self' = 'Ideology',
                                           'factor(employment)Blue collar' = 'Blue collar',
                                           'factor(employment)White collar' = 'White collar',
                                           'subjective_class' = 'Subjective Class',
                                           'living_area' = 'Rural / urban',
                                           'local_interest' = 'Interest in local politics',
                                           'national_interest' = 'Interest in national politics',
                                           'european_interest' = 'Interest in European politics'),
                              title = 'Effect of the UK vaccine rollout on EU attitudes',
                              notes = 'Note: Robust standard errors in parentheses. Reference category for employment is "Unemployment". All models include country fixed effects.',
                              output = '~/Replication materials/Tables/app_table_5.png')

#### Table 6: Effect of the UK vaccine rollout on attitudes towards European decision-making ####
## Model building for DVs: EU health issues, social security, education, equality, climate, jobs, digitalisation, and working conditions
model_eu_health_issues <- lm_robust(eu_health_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                              data = data,
                                              fixed_effects = country,
                                              se_type = "stata")
model_socialsecurity <- lm_robust(eu_ss_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                            data = data,
                                            fixed_effects = country,
                                            se_type = "stata")
model_education <- lm_robust(eu_education_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                       data = data,
                                       fixed_effects = country,
                                       se_type = "stata")
model_equality <- lm_robust(eu_equality_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                      data = data,
                                      fixed_effects = country,
                                      se_type = "stata")
model_climate <- lm_robust(eu_climate_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                     data = data,
                                     fixed_effects = country,
                                     se_type = "stata")
model_jobs <- lm_robust(eu_jobs_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                  data = data,
                                  fixed_effects = country,
                                  se_type = "stata")
model_digitalisation <- lm_robust(eu_digitalisation_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                            data = data,
                                            fixed_effects = country,
                                            se_type = "stata")
model_workingconditions <- lm_robust(eu_workingconditions_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                               data = data,
                                               fixed_effects = country,
                                               se_type = "stata")

## Making a list of all Specific variables models
models_specific <- list("Health" = model_eu_health_issues,
                                  "Social security" = model_socialsecurity,
                                  "Education" = model_education,
                                  "Equality" = model_equality,
                                  "Climate" = model_climate,
                                  "Working conditions" = model_workingconditions,
                                  "Jobs" = model_jobs,
                                  "Digitalisation" = model_digitalisation)

## Building the regression table of specific variables models and saving Table 6 as 'app_table_6.png'
t_models_specific <- modelsummary(models_specific,
                                  stars = TRUE,
                                  coef_map = c('treat2' = 'Treatment',
                                               'gender' = 'Gender',
                                               'age' = 'Age',
                                               'eduy' = 'Education',
                                               'lr_self' = 'Ideology',
                                               'factor(employment)Blue collar' = 'Blue collar',
                                               'factor(employment)White collar' = 'White collar',
                                               'subjective_class' = 'Subjective Class',
                                               'living_area' = 'Rural / urban',
                                               'local_interest' = 'Interest in local politics',
                                               'national_interest' = 'Interest in national politics',
                                               'european_interest' = 'Interest in European politics'),
                                  title = 'Effect of the UK vaccine rollout on EU attitudes',
                                  notes = 'Note: Robust standard errors in parentheses. Reference category for employment is "Unemployment". All models include country fixed effects.',
                                  output = '~/Replication materials/Tables/app_table_6.png')


#### Table 7: Effect of the UK vaccine rollout on EU attitudes - moderating effect of regions ####

## Model building for interaction with region with DVs: EU support, EU benefits, EU positive image, EU right direction, EU democracy, EU health issues, EU coordination, and EU future coordination
model_eu_support_regions <- lm_robust(eu_support ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                      data = data,
                                      se_type = "stata")

model_eu_beneficial_regions <- lm_robust(eu_benefits ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                         data = data,
                                         se_type = "stata")

model_eu_positive_regions <- lm_robust(eu_positive ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                       data = data,
                                       se_type = "stata")

model_eu_health_issues_regions <- lm_robust(eu_health_issues ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                            data = data,
                                            se_type = "stata")

model_eu_direction_regions <- lm_robust(eu_right_direction ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                        data = data,
                                        se_type = "stata")

model_eu_democracy_regions <- lm_robust(eu_democracy ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                        data = data,
                                        se_type = "stata")

model_eu_coordination_regions <- lm_robust(eu_help ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                           data = data,
                                           se_type = "stata")

model_eu_future_coordination_regions <- lm_robust(eu_coord ~ treat2*factor(europe) + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                  data = data,
                                                  se_type = "stata")

## Making a list of all Region interaction models
models_regions <- list("EU support" = model_eu_support_regions,
                       "EU beneficial" = model_eu_beneficial_regions,
                       "Positive image" = model_eu_positive_regions,
                       "EU right direction" = model_eu_direction_regions,
                       "EU democracy" = model_eu_democracy_regions,
                       "Health decisions" = model_eu_health_issues_regions,
                       "EU crisis coordination" = model_eu_coordination_regions,
                       "EU future crisis coordination" = model_eu_future_coordination_regions)

## ## Building the regression table of specific variables models and saving Table 7 as 'app_table_7.png'
t_models_regions <- modelsummary(models_regions,
                                 stars = TRUE,
                                 coef_map = c('treat2:factor(europe)eastern' = 'Treatment x Eastern region',
                                              'treat2:factor(europe)western' = 'Treatment x Western region',
                                              'treat2:factor(europe)northern' = 'Treatment x Northern region',
                                              'treat2:factor(europe)southern' = 'Treatment x Southern region',
                                              '(Intercept)' = 'Intercept'),
                                 notes = 'Robust standard errors in parentheses. Controls include: age, gender, education, left-right self-placement, emmployment status, subjective class, living area and interest in local, national and European politics.',
                                 title = 'Effect of the UK vaccine rollout on EU attitudes - moderating effect of regions',
                                 output = "~/Replication materials/Tables/app_table_7.png")

#### Figure 6: Predicted values of EU attitudes by region of the respondent ####
## To be able to plot predicted values of models interacting treatment with region, we replicate the models using a simple lm() regression.
model_eu_benefits_lm_regions <- lm(eu_benefits ~ treat2*europe + age + gender + eduy + lr_self + employment + subjective_class + living_area + local_interest + national_interest + european_interest,
                                   data = data)
model_eu_health_issues_lm_regions <- lm(eu_health_issues ~ treat2*europe + age + gender + eduy + lr_self + employment + subjective_class + living_area + local_interest + national_interest + european_interest,
                                        data = data)
## Building plot for Figure 6A
p_lm_eu_benefits_regions <- plot_predictions(model_eu_benefits_lm_regions,
                                             condition = c("treat2", "europe")) +
  labs(title = "",
       x = "Treatment status",
       y = "Predicted values",
       color = "Region in \nEurope") +
  scale_x_discrete(labels = c("Before vaccination news", "After vaccination news")) +
  scale_colour_discrete(labels = c("Central", "Eastern", "Northern", "Southern", "Western")) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 20))

## Saving Figure 6A as 'app_figure_6a.png'
ggsave("~/Replication materials/Figures/app_figure_6a.png", 
       plot = p_lm_eu_benefits_regions, 
       dpi = 600, 
       width = 8, 
       height = 7)


## Building plot for Figure 6B
p_lm_eu_health_issues_regions <- plot_predictions(model_eu_health_issues_lm_regions,
                                                  condition = c("treat2", "europe")) +
  labs(title = "",
       x = "Treatment status",
       y = "Predicted values") +
  scale_x_discrete(labels = c("Before vaccination news", "After vaccination news")) +
  scale_colour_discrete(labels = c("Central", "Eastern", "Northern", "Southern", "Western")) +
  theme_classic() +
  theme(axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        legend.position="none")

## Saving Figure 6B as 'app_figure_6b.png'
ggsave("~/Replication materials/Figures/app_figure_6b.png", 
       plot = p_lm_eu_health_issues_regions, 
       dpi = 600, 
       width = 8, 
       height = 7)

#### Table 8: Effect of the UK vaccine rollout on EU attitudes - interaction with polarisation ####

## Model building for interaction with polarisation with DVs: EU support, EU benefits, EU health issues, EU coordination
model_eu_support_polarisation <- lm_robust(eu_support ~ treat2*elite_ideo + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                   data = data,
                                                   se_type = "stata")
model_eu_beneficial_polarisation <- lm_robust(eu_benefits ~ treat2*elite_ideo + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                      data = data,
                                                      se_type = "stata")

model_eu_health_issues_polarisation <- lm_robust(eu_health_issues ~ treat2*elite_ideo + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                         data = data,
                                                         se_type = "stata")

model_eu_coordination_polarisation <- lm_robust(eu_help ~ treat2*elite_ideo + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                        data = data,
                                                        se_type = "stata")

## Making a list of all polarisation models
models_polarisation <- list("EU support" = model_eu_support_polarisation,
                            "EU benefits" = model_eu_beneficial_polarisation,
                            "EU health decisions" = model_eu_health_issues_polarisation,
                            "EU crisis coordination" = model_eu_coordination_polarisation)

## Building regression table for Table 8 and saving it as 'app_table_8.png'
t_models_polarisation <- modelsummary(models_polarisation, 
                                      stars = TRUE,
                                      coef_map = c(
                                        'treat2:elite_ideo' = 'Treatment x Elite polarisation',
                                        'treat2' = 'Treatment',
                                        'elite_ideo' = 'Elite polarisation',
                                        '(Intercept)' = 'Intercept'),
                                      title = 'Effect of the UK vaccine rollout on EU attitudes - interaction with polarisation',
                                      notes = 'Note: Robust standard errors in parentheses',
                                      output = '~/Replication materials/Tables/app_table_8.png')


#### Table 9: Effect of the UK vaccine rollout on EU attitudes - moderating effect of hospital occupancy rates (normalised between 0 and 1) ####
## Model building for interaction with polarisation with DVs: EU support, EU benefits, EU positive image, EU right direction, EU democracy, EU health issues, EU coordination, and EU future coordination
model_eu_support_hospitals <- lm_robust(eu_support ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                  data = data,
                                                  fixed_effects = country,
                                                  se_type = "stata")

model_eu_beneficial_hospitals <- lm_robust(eu_benefits ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                     data = data,
                                                     fixed_effects = country,
                                                     se_type = "stata")

model_eu_positive_hospitals <- lm_robust(eu_positive ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                   data = data,
                                                   fixed_effects = country,
                                                   se_type = "stata")

model_eu_health_issues_hospitals <- lm_robust(eu_health_issues ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                        data = data,
                                                        fixed_effects = country,
                                                        se_type = "stata")

model_eu_right_direction_hospitals <- lm_robust(eu_right_direction ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                          data = data,
                                                          fixed_effects = country,
                                                          se_type = "stata")

model_eu_democracy_hospitals <- lm_robust(eu_democracy ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                    data = data,
                                                    fixed_effects = country,
                                                    se_type = "stata")

model_eu_coordination_hospitals <- lm_robust(eu_help ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                       data = data,
                                                       fixed_effects = country,
                                                       se_type = "stata")

model_eu_future_coordination_hospitals <- lm_robust(eu_coord ~ treat2*hospital_occupancy + age + gender + eduy + lr_self + factor(employment) + subjective_class + living_area + local_interest + national_interest + european_interest,
                                                              data = data,
                                                              fixed_effects = country,
                                                              se_type = "stata")

## Making a list of all models interacting with COVID-19-related hospitalisations
models_hospitals <- list("EU support" = model_eu_support_hospitals,
                         "EU beneficial" = model_eu_beneficial_hospitals,
                         "Positive image" = model_eu_positive_hospitals,
                         "EU right direction" = model_eu_right_direction_hospitals,
                         "EU democracy" = model_eu_democracy_hospitals,
                         "Health decisions" = model_eu_health_issues_hospitals,
                         "EU crisis coordination" = model_eu_coordination_hospitals,
                         "EU future crisis coordination" = model_eu_future_coordination_hospitals)

## Building regression for Table 9 and saving it as 'app_table_9.png'
t_models_hospitals <- modelsummary(models_hospitals,
                                   stars = TRUE,
                                   coef_map = c('treat2:hospital_occupancy' = 'Treatment x hospital occupancy rates',
                                                'hospital_occupancy' = 'Hospital occupancy rates',
                                                'treat2' = 'Treatment', 
                                                '(Intercept)' = 'Intercept'),
                                   title = 'Effect of the UK vaccine rollout on EU attitudes - moderating effect of hospital occupancy rates (normalised between 0 and 1) - portugal + FE',
                                   notes = c('Note: Robust standard errors in parentheses. Controls include: age, gender, education, left-right self-placement, employment status, subjective class, living area and interest in local, national and European politics.'),
                                   output = "~/Replication materials/Tables/app_table_9.png")


#### Table 10: Models Excluding Different Countries on Health Decisions, EU Coordination, and Future EU Coordination ####

# List of countries to exclude
exclude_list <- list(
  "Austria", 
  "Cyprus", 
  "Netherlands", 
  "Spain",
  c("Austria", "Cyprus"),
  c("Austria", "Netherlands"),
  c("Austria", "Spain"),
  c("Cyprus", "Netherlands"),
  c("Cyprus", "Spain"),
  c("Netherlands", "Spain"),
  c("Austria", "Cyprus", "Netherlands"),
  c("Austria", "Cyprus", "Spain"),
  c("Cyprus", "Netherlands", "Spain"),
  c("Austria", "Netherlands", "Spain"),
  c("Austria", "Cyprus", "Netherlands", "Spain")
)

# Function to run models excluding specific countries
run_models <- function(exclude) {
  data_filtered <- data |> filter(!country %in% exclude)
  
  model_health <- lm_robust(eu_health_issues ~ treat2 + age + gender + eduy + lr_self + factor(employment) + living_area + national_interest + european_interest + local_interest,
                            data = data_filtered, fixed_effects = country, se_type = "stata")
  
  model_coordination <- lm_robust(eu_help ~ treat2 + age + gender + eduy + lr_self + factor(employment) + living_area + national_interest + european_interest + local_interest,
                                  data = data_filtered, fixed_effects = country, se_type = "stata")
  
  model_future <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + lr_self + factor(employment) + living_area + national_interest + european_interest + local_interest,
                            data = data_filtered, fixed_effects = country, se_type = "stata")
  
  return(list("health" = model_health, "coordination" = model_coordination, "future" = model_future))
}

# Applying the function to all country exclusions
model_results <- lapply(exclude_list, run_models)
# Extracting coefficient and standard error
extract_coeff_se <- function(model) {
  coef <- coef(model)["treat2"]
  se <- sqrt(vcov(model)["treat2", "treat2"])
  return(c(coef, se))
}

# Creating a data frame to hold the results
results_df <- data.frame(Country_Exclusion = character(), 
                         Health_Decisions_Coef = numeric(), Health_Decisions_SE = numeric(),
                         EU_Coordination_Coef = numeric(), EU_Coordination_SE = numeric(),
                         Future_Coordination_Coef = numeric(), Future_Coordination_SE = numeric(),
                         stringsAsFactors = FALSE)

# Adding row names for the exclusions
country_labels <- c("Excluding Austria", "Excluding Cyprus", "Excluding Netherlands", 
                    "Excluding Spain", "Excluding Austria and Cyprus", 
                    "Excluding Austria and Netherlands", "Excluding Austria and Spain", 
                    "Excluding Cyprus and Netherlands", "Excluding Cyprus and Spain", 
                    "Excluding Netherlands and Spain", "Excluding Austria, Cyprus, Netherlands",
                    "Excluding Austria, Cyprus, Spain", "Excluding Cyprus, Netherlands, Spain",
                    "Excluding Austria, Netherlands, Spain", "Excluding all four")

# Loop through model results and extract coefficients & SEs
for (i in seq_along(model_results)) {
  model_set <- model_results[[i]]
  
  health_results <- extract_coeff_se(model_set$health)
  coord_results <- extract_coeff_se(model_set$coordination)
  future_results <- extract_coeff_se(model_set$future)
  
  # Add row to the results data frame
  results_df <- rbind(results_df, data.frame(
    Country_Exclusion = country_labels[i],
    Health_Decisions_Coef = health_results[1], Health_Decisions_SE = health_results[2],
    EU_Coordination_Coef = coord_results[1], EU_Coordination_SE = coord_results[2],
    Future_Coordination_Coef = future_results[1], Future_Coordination_SE = future_results[2]
  ))
}

# Round the coefficients and SEs for readability
results_df <- results_df %>%
  mutate(across(contains("Coef") | contains("SE"), round, 3))


# Adding significance stars based on p-values
add_significance <- function(coef, se) {
  p_value <- 2 * (1 - pnorm(abs(coef / se)))
  if (p_value < 0.001) return("***")
  if (p_value < 0.01) return("**")
  if (p_value < 0.05) return("*")
  if (p_value < 0.1) return("+")
  return("")
}

# Adding significance stars to the table using mapply (vectorized operation)
results_df <- results_df %>%
  mutate(Health_Decisions = paste0(Health_Decisions_Coef, 
                                   mapply(add_significance, Health_Decisions_Coef, Health_Decisions_SE),
                                   "\n(", Health_Decisions_SE, ")"),
         EU_Coordination = paste0(EU_Coordination_Coef, 
                                  mapply(add_significance, EU_Coordination_Coef, EU_Coordination_SE),
                                  "\n(", EU_Coordination_SE, ")"),
         Future_Coordination = paste0(Future_Coordination_Coef, 
                                      mapply(add_significance, Future_Coordination_Coef, Future_Coordination_SE),
                                      "\n(", Future_Coordination_SE, ")"))

# Selecting and reorder the columns for output
results_table <- results_df %>%
  select(Country_Exclusion, Health_Decisions, EU_Coordination, Future_Coordination)


# Creating the table using kableExtra
table_html <- results_table %>%
  kbl(caption = "Models Excluding Different Countries on Health Decisions, EU Coordination, and Future EU Coordination") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Save the Table 10 as HTML as 'app_table_10.html'
save_kable(table_html, "Tables/app_table_10.html")


#### Table 11: Results of Different Matching Strategies on Health Decisions, EU Coordination, and Future EU Coordination ####

## Matching 1 - Matching by Age, gender, education and living area

# Filtering data by variables of interest and deleting missing observations
dat_matching1 <- data |>
  select(c("treat2", "age", "gender", "eduy", "living_area", "country",
           "eu_support", "eu_positive", "eu_benefits", "eu_health_issues", "eu_right_direction", "eu_democracy", "eu_help", "eu_coord", "europe"))  |> 
  na.omit()

# Matching observations by treatment assignment by covariates
model_matching1 <- matchit(treat2 ~ age + gender + eduy + living_area,
                           data = dat_matching1)

matched_data_matching1 <- match.data(model_matching1)

matching1_health <- lm_robust(eu_health_issues~ treat2 + age + gender + eduy + living_area,
                              data = matched_data_matching1,
                              fixed_effects = country,
                              se_type = "stata")

matching1_coordination <- lm_robust(eu_help ~ treat2 + age + gender + eduy + living_area,
                                    data = matched_data_matching1,
                                    fixed_effects = country,
                                    se_type = "stata")

matching1_future <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + living_area,
                              data = matched_data_matching1,
                              fixed_effects = country,
                              se_type = "stata")


# Extract p-values, standard errors and coefficients
p_matching1_health <- summary(matching1_health)$coefficients["treat2", "Pr(>|t|)"]
se_matching1_health <- summary(matching1_health)$coefficients["treat2", "Std. Error"]
coef_matching1_health <- coef(matching1_health)["treat2"]

p_matching1_coordination <- summary(matching1_coordination)$coefficients["treat2", "Pr(>|t|)"]
se_matching1_coordination <- summary(matching1_coordination)$coefficients["treat2", "Std. Error"]
coef_matching1_coordination <- coef(matching1_coordination)["treat2"]

p_matching1_future <- summary(matching1_future)$coefficients["treat2", "Pr(>|t|)"]
se_matching1_future <- summary(matching1_future)$coefficients["treat2", "Std. Error"]
coef_matching1_future <- coef(matching1_future)["treat2"]




## Matching 2 - Adding employment status
dat_matching2 <- data |>
  select(c("treat2", "age", "gender", "eduy", "employment_r", "living_area", "country",
           "eu_support", "eu_positive", "eu_benefits", "eu_health_issues", "eu_right_direction", "eu_democracy", "eu_help", "eu_coord", "europe")) |>
  na.omit()


model_matching2 <- matchit(treat2 ~ age + gender + eduy + living_area + employment_r,
                           data = dat_matching2)

matched_data_matching2 <- match.data(model_matching2)


matching2_health <- lm_robust(eu_health_issues~ treat2 + age + gender + eduy + employment_r + living_area,
                              data = matched_data_matching2,
                              fixed_effects = country,
                              se_type = "stata")
matching2_coordination <- lm_robust(eu_help ~ treat2 + age + gender + eduy + employment_r + living_area,
                                    data = matched_data_matching2,
                                    fixed_effects = country,
                                    se_type = "stata")

matching2_future <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + employment_r + living_area,
                              data = matched_data_matching2,
                              fixed_effects = country,
                              se_type = "stata")

# Extract p-values, standard errors and coefficients
p_matching2_health <- summary(matching2_health)$coefficients["treat2", "Pr(>|t|)"]
se_matching2_health <- summary(matching2_health)$coefficients["treat2", "Std. Error"]
coef_matching2_health <- coef(matching2_health)["treat2"]

p_matching2_coordination <- summary(matching2_coordination)$coefficients["treat2", "Pr(>|t|)"]
se_matching2_coordination <- summary(matching2_coordination)$coefficients["treat2", "Std. Error"]
coef_matching2_coordination <- coef(matching2_coordination)["treat2"]

p_matching2_future <- summary(matching2_future)$coefficients["treat2", "Pr(>|t|)"]
se_matching2_future <- summary(matching2_future)$coefficients["treat2", "Std. Error"]
coef_matching2_future <- coef(matching2_future)["treat2"]



## Matching 3 - Adding matching by left-right self-placement
dat_matching3 <- data |>
  select(c("treat2", "age", "gender", "eduy", "employment_r", "lr_self", "living_area", "country",
           "eu_support", "eu_positive", "eu_benefits", "eu_health_issues", "eu_right_direction", "eu_democracy", "eu_help", "eu_coord", "europe")) |>
  na.omit()

model_matching3 <- matchit(treat2 ~ age + gender + eduy + living_area + employment_r + lr_self,
                           data = dat_matching3)

matched_data_matching3 <- match.data(model_matching3)


matching3_health <- lm_robust(eu_health_issues~ treat2 + age + gender + eduy + employment_r + living_area + lr_self,
                              data = matched_data_matching3,
                              fixed_effects = country,
                              se_type = "stata")
matching3_coordination <- lm_robust(eu_help ~ treat2 + age + gender + eduy + employment_r + living_area + lr_self,
                                    data = matched_data_matching3,
                                    fixed_effects = country,
                                    se_type = "stata")

matching3_future <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + employment_r + living_area + lr_self,
                              data = matched_data_matching3,
                              fixed_effects = country,
                              se_type = "stata")


# Extract p-values, standard errors and coefficients
p_matching3_health <- summary(matching3_health)$coefficients["treat2", "Pr(>|t|)"]
se_matching3_health <- summary(matching3_health)$coefficients["treat2", "Std. Error"]
coef_matching3_health <- coef(matching3_health)["treat2"]

p_matching3_coordination <- summary(matching3_coordination)$coefficients["treat2", "Pr(>|t|)"]
se_matching3_coordination <- summary(matching3_coordination)$coefficients["treat2", "Std. Error"]
coef_matching3_coordination <- coef(matching3_coordination)["treat2"]

p_matching3_future <- summary(matching3_future)$coefficients["treat2", "Pr(>|t|)"]
se_matching3_future <- summary(matching3_future)$coefficients["treat2", "Std. Error"]
coef_matching3_future <- coef(matching3_future)["treat2"]


## Matching 4 - Adding matching by country

dat_matching4 <- data |>
  select(c("treat2", "age", "gender", "eduy", "lr_self", "employment_r", "country", "living_area",
           "eu_support", "eu_positive", "eu_benefits", "eu_health_issues", "eu_right_direction", "eu_democracy", "eu_help", "eu_coord", "europe")) |>
  na.omit()



model_matching4 <- matchit(treat2 ~ age + gender + eduy + living_area + employment_r + lr_self + country,
                           data = dat_matching4)

matched_data_matching4 <- match.data(model_matching4)


matching4_health <- lm_robust(eu_health_issues~ treat2 + age + gender + eduy + employment_r + living_area + lr_self + country,
                              data = matched_data_matching4,
                              se_type = "stata")
matching4_coordination <- lm_robust(eu_help ~ treat2 + age + gender + eduy + employment_r + living_area + lr_self + country,
                                    data = matched_data_matching4,
                                    se_type = "stata")

matching4_future <- lm_robust(eu_coord ~ treat2 + age + gender + eduy + employment_r + living_area + lr_self + country,
                              data = matched_data_matching4,
                              se_type = "stata")


# Extract p-values, standard errors and coefficients
p_matching4_health <- summary(matching4_health)$coefficients["treat2", "Pr(>|t|)"]
se_matching4_health <- summary(matching4_health)$coefficients["treat2", "Std. Error"]
coef_matching4_health <- coef(matching4_health)["treat2"]

p_matching4_coordination <- summary(matching4_coordination)$coefficients["treat2", "Pr(>|t|)"]
se_matching4_coordination <- summary(matching4_coordination)$coefficients["treat2", "Std. Error"]
coef_matching4_coordination <- coef(matching4_coordination)["treat2"]

p_matching4_future <- summary(matching4_future)$coefficients["treat2", "Pr(>|t|)"]
se_matching4_future <- summary(matching4_future)$coefficients["treat2", "Std. Error"]
coef_matching4_future <- coef(matching4_future)["treat2"]

## Adding all data to a data frame to form Table 11
results_matching <- data.frame(
  Matching = c("Matching 1", "Matching 2", "Matching 3", "Matching 4"),
  Health_Decisions = c(
    paste0(round(coef_matching1_health, 3), ifelse(p_matching1_health < 0.001, "***", 
                                                   ifelse(p_matching1_health < 0.01, "**", 
                                                          ifelse(p_matching1_health < 0.05, "*", 
                                                                 ifelse(p_matching1_health < 0.1, "+", "")))), " (", round(se_matching1_health, 3), ")", ""),
    paste0(round(coef_matching2_health, 3), ifelse(p_matching2_health < 0.001, "***", 
                                                   ifelse(p_matching2_health < 0.01, "**", 
                                                          ifelse(p_matching2_health < 0.05, "*", 
                                                                 ifelse(p_matching2_health < 0.1, "+", "")))), " (", round(se_matching2_health, 3), ")", ""),
    paste0(round(coef_matching3_health, 3), ifelse(p_matching3_health < 0.001, "***", 
                                                   ifelse(p_matching3_health < 0.01, "**", 
                                                          ifelse(p_matching3_health < 0.05, "*", 
                                                                 ifelse(p_matching3_health < 0.1, "+", "")))), " (", round(se_matching3_health, 3), ")", ""),
    paste0(round(coef_matching4_health, 3), ifelse(p_matching4_health < 0.001, "***", 
                                                   ifelse(p_matching4_health < 0.01, "**", 
                                                          ifelse(p_matching4_health < 0.05, "*", 
                                                                 ifelse(p_matching4_health < 0.1, "+", "")))), " (", round(se_matching4_health, 3), ")", "")
  ),
  EU_Coordination = c(
    paste0(round(coef_matching1_coordination, 3), ifelse(p_matching1_coordination < 0.001, "***", 
                                                         ifelse(p_matching1_coordination < 0.01, "**", 
                                                                ifelse(p_matching1_coordination < 0.05, "*", 
                                                                       ifelse(p_matching1_coordination < 0.1, "+", "")))), " (", round(se_matching1_coordination, 3), ")", ""),
    paste0(round(coef_matching2_coordination, 3), ifelse(p_matching2_coordination < 0.001, "***", 
                                                         ifelse(p_matching2_coordination < 0.01, "**", 
                                                                ifelse(p_matching2_coordination < 0.05, "*", 
                                                                       ifelse(p_matching2_coordination < 0.1, "+", "")))), " (", round(se_matching2_coordination, 3), ")", ""),
    paste0(round(coef_matching3_coordination, 3), ifelse(p_matching3_coordination < 0.001, "***", 
                                                         ifelse(p_matching3_coordination < 0.01, "**", 
                                                                ifelse(p_matching3_coordination < 0.05, "*", 
                                                                       ifelse(p_matching3_coordination < 0.1, "+", "")))), " (", round(se_matching3_coordination, 3), ")", ""),
    paste0(round(coef_matching4_coordination, 3), ifelse(p_matching4_coordination < 0.001, "***", 
                                                         ifelse(p_matching4_coordination < 0.01, "**", 
                                                                ifelse(p_matching4_coordination < 0.05, "*", 
                                                                       ifelse(p_matching4_coordination < 0.1, "+", "")))), " (", round(se_matching4_coordination, 3), ")", "")
  ),
  EU_Future_Coordination = c(
    paste0(round(coef_matching1_future, 3), ifelse(p_matching1_future < 0.001, "***", 
                                                   ifelse(p_matching1_future < 0.01, "**", 
                                                          ifelse(p_matching1_future < 0.05, "*", 
                                                                 ifelse(p_matching1_future < 0.1, "+", "")))), " (", round(se_matching1_future, 3), ")", ""),
    paste0(round(coef_matching2_future, 3), ifelse(p_matching2_future < 0.001, "***", 
                                                   ifelse(p_matching2_future < 0.01, "**", 
                                                          ifelse(p_matching2_future < 0.05, "*", 
                                                                 ifelse(p_matching2_future < 0.1, "+", "")))), " (", round(se_matching2_future, 3), ")", ""),
    paste0(round(coef_matching3_future, 3), ifelse(p_matching3_future < 0.001, "***", 
                                                   ifelse(p_matching3_future < 0.01, "**", 
                                                          ifelse(p_matching3_future < 0.05, "*", 
                                                                 ifelse(p_matching3_future < 0.1, "+", "")))), " (", round(se_matching3_future, 3), ")", ""),
    paste0(round(coef_matching4_future, 3), ifelse(p_matching4_future < 0.001, "***", 
                                                   ifelse(p_matching4_future < 0.01, "**", 
                                                          ifelse(p_matching4_future < 0.05, "*", 
                                                                 ifelse(p_matching4_future < 0.1, "+", "")))), " (", round(se_matching4_future, 3), ")", "")
  )
)


# Creating and styling the table
table_html <- results_matching %>%
  kbl(caption = "Results of Different Matching Strategies on Health Decisions, EU Coordination, and Future EU Coordination") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed")) %>%
  add_header_above(c(" " = 1, "Health Decisions" = 1, "EU Coordination" = 1, "EU Future Coordination" = 1)) %>%
  column_spec(4, bold = TRUE, italic = TRUE) # To style the description column

# Save the table as HTML as 'app_table_11.html'
save_kable(table_html, "Tables/app_table_11.html")



#### Figure 7: Factors found on EU specific attitudes ####

set.seed(123)

## Filtering data by variables of interest and deleting all observations with missing values
dat_fa_2 <- data %>%
  mutate(eu_health_issues = as.numeric(eu_health_issues),
         eu_ss_issues = as.numeric(eu_ss_issues),
         eu_education_issues = as.numeric(eu_education_issues),
         eu_equality_issues = as.numeric(eu_equality_issues),
         eu_climate_issues = as.numeric(eu_climate_issues),
         eu_jobs_issues = as.numeric(eu_jobs_issues),
         eu_digitalisation_issues = as.numeric(eu_digitalisation_issues),
         eu_workingconditions_issues = as.numeric(eu_workingconditions_issues),
         eu_benefits = as.numeric(eu_benefits)) %>%
  select(c(eu_health_issues, eu_ss_issues, eu_education_issues, eu_equality_issues,
           eu_climate_issues, eu_jobs_issues, eu_digitalisation_issues, eu_workingconditions_issues, eu_benefits)) %>%
  na.omit() 

## Principal Axis factor analysis 
fa.none <- fa(r = dat_fa_2,
              nfactors = 3, # three factors
              fm = "pa", # principal axis
              max.iter = 100, # iterations
              SMC = FALSE)

## Visual representation of the three factors and the elements that compose them
fa.diagram(fa.none)

## Saving output of previous step as 'app_figure_7.png'
dev.print(device = png,
          filename = "~/Replication materials/Figures/app_figure_7.png",
          width = 600,              
          height = 550              
)
